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A general stability criterion is derived for the ground states of the Gross-Pitaevskii equation, which 
describes attractive Bose-Einstein condensates confined in a magnetic trap. These ground states are 
shown to avoid the collapse in finite time and are proven to be stable in two and three spatial dimensions. 
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I. INTRODUCTION 

Experimental observation of Bose-Einstein condensa- 
tion (BEC) in ultracold atomic clouds has stimulated 
a new direction in the study of macroscopic quantum phe- 
nomena. Basically, the interaction between two confined 
bosons in a condensate is determined by the s-wave scat- 
tering length clq and it can be either repulsive (ao > 0) 
or attractive (ao < 0). Although first BEC experiments 
were commonly realized with gases promoting a posi- 
tive scattering length, trapped 7 Li atom gases, which 
are characterized by a negative scattering length, have 
raised an increasing interest justified by the rich and 
complex dynamics mixing instability and generation of 
solitonlike structures, which substantially alters the for- 
mation of condensates. Furthermore, novel experimental 
results {I suggested the possibility of using so-called Fes- 
hbach resonances to continuously detune ao from positive 
to negative values by means of an external magnetic field, 
which brings a new insight into the experimental realiza- 
tion of BECs with attractive interactions. 

The dynamical behaviors of gases with negative scat- 
tering length are especially interesting, because in the 
absence of trapping the condensate is described by the 
nonlinear Schrodinger (NLS) equation, whose localized, 
multi-dimensional solutions are unstable and may col- 
lapse in finite time (for a review see, e.g., Ref. (J]). Sim- 
ilarly, collapse also occurs in confined condensates with 
attractive interactions, as emphasized by many theoret- 
ical works HJ-Q], and sequences of collapse events have 
been experimentally detected in BECs of 7 Li gas ||, in 
which the condensate was observed to shrink on time 
scales of the trap oscillation period. General conditions 
for collapse follow from the virial theorem applied to the 
Gross-Pitaevskii (GP) equation 

ih^ = +^f(r 2 + 7 V)V> + WIV, (1) 

where ip denotes the BEC wavefunction, V 2 is the three- 
dimensional Laplacian, 7 measures the anisotropy ratio 



of the condensate and Uq = Airh 2 ao j 'm characterizes the 
two-body interaction (771 is the atom mass and u) the 
trap frequency). These conditions imply the existence 
of a critical value, N CT , that the total number of atoms 
N = J I tp 1 2 dr must exceed for initiating the collapse. For 
N below critical, attractive condensates are now known 
to be capable of relaxing to stationary states of the GP 
equation, reading as ip{r,t) — 4>{r) exp (—i/it/K), where 
/i has the meaning of a chemical potential. These solu- 
tions consist of solitonlike ground states, whose stability 
was recently analyzed by means of numerical continua- 
tion methods |9| . ft was in particular revealed that as N 
approaches N CT from below, the condensate modeled by 
this ground state becomes more unstable and no station- 
ary condensate forms for N > N cr . Besides, a new inves- 
tigation |l(| showed that BEC with large negative scat- 
tering length can be modulationally unstable and decay 
into periodic fringe patterns of parallel solitons, which 
again stresses the important role of the latter in this field. 

In spite of these recent results, some fundamental 
points, however, remain unsolved. First, conditions for 
the existence of stationary condensates have never been 
established for high-dimensional systems and we still do 
not know how they may be compatible with the condi- 
tions for collapse, which mainly follow from the Hamil- 
tonian properties of the GP equation. Second, we ig- 
nore how trapped condensates can evolve towards ground 
states. Third, for practical uses, it is highly desirable 
to know whether stability of stationary condensates can 
be inferred from their number of particles, and from 
this number only. The aim of this paper is to answer 
these questions and to clear up the link between the 
non-occurrence of collapse with respect to the stability 
of localized, stationary ground states when ao < 0. In 
the following, we review the conditions required for the 
collapse of two-dimensional (2D) and three-dimensional 
(3D) condensates and present a rigorous derivation of 
the stability criterion for the ground-state solutions of 
the GP equation. 3D isotropic BECs are analyzed in 



the limit 7 



1, for which radially-symmetric solu- 



tions only depend on the radius r = \J x 1 + y 2 + z 2 in 
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Eq. (@). In 2D, the model (@) can be employed to 
describe the dynamics of anisotropic condensates hav- 
ing a pancake-shape geometry with j 2 3> 1. In this 
configuration, ip = ip'(x,y,t)<f>(z) exhibits a longitu- 
dinal component, (f)(z), frozen on the Gaussian shape 
(j)(z) = ( 7 / 7 r) 1 / 4 exp(-7z 2 /2) The condensate 

dynamics is then basically governed by the 2D version 
of Eq. (jl]) with radius r = \J x 2 + y 2 . ID BECs are 
disregarded, since related results are already available in 
the current literature devoted to optical solitons |E| . In- 
stead, we focus our attention on situations allowing for 
collapse, which basically concerns high dimension num- 
bers fl3[| . As collapse properties depend on more than 
the dimensionality rather than the geometry of nonlinear 
objects, we only address the stability of isotropic conden- 
sates. 



H = X - Y + ^N(r 2 ), (4) 

where 

X ee J \Vu\ 2 dr, y=\j (5) 

fulfils conditions making the right-hand side of ([?]) nega- 
tive. 

Two analytical estimates are usually employed to de- 
termine necessary requirements on the particle number 
for initiating the collapse |p^ , p^| , namely, the "uncer- 
tainty" inequality 

N<(2/D) 2 X(r 2 ), (6) 



II. COLLAPSE VS. GROUND STATES 

^From the viewpoint of the nonlinear dynamics of lo- 
calized modes, studying multi-dimensional trapped con- 
densates is equivalent to investigating the solutions of the 
Z)-dimensional equation 

id t u + V 2 u- -^-r 2 u + \u\ 2 u = 0, (2) 

where dt is a partial derivative in time and V 2 = 
r 1 ~ D d r r D ~ 1 d r . In Eq. @, the time and space vari- 
ables are measured in units {h/muj) 1 / 2 and 2/u, respec- 
tively; the new wavefunction u scales as (Huj/2Uo)^ 1 ^ 2 , 
and f2o — 2 stands for the trap frequency u>, which is for- 
mally kept undefined for technical convenience. Briefly 
speaking, collapse occurs as the mean-square radius of 
the wavefunction, (r 2 ) = (1/N) J r 2 \u\ 2 dr, where N — 
J \u\ 2 dr accounts for the rescaled particle number, tends 
to zero in finite time, while the amplitude of u(r,t) di- 
verges in the same limit. To be triggered off, this singu- 
lar phenomenon requires special conditions on the ini- 
tial data R 2 = (r 2 )| t=0 and R 2 . = d t (r 2 )\ t =o, among 
which we henceforth consider wavefunctions with no ini- 
tial divergence and set R = 0, for the sake of clarity. 
These requirements are basically inferred from the so- 
called "virial" relation 

9tV) = ^ + / \u\'dr - 4fi 2 (r 2 ), (3) 

which is merely derived by multiplying Eq.(j2|) with (r 2 u*) 
and (r • Vu*) (* means complex conjugate) and by com- 
bining the imaginary and real parts of the results, respec- 
tively d-QJlH . Equation (||) then shows that the mean- 
square radius (r 2 (t)) inevitably vanishes at a finite time, 
t c , whenever D > 2 and when the conserved Hamiltonian 
integral 



from which the gradient norm X blows up as (r 2 ) — ► 0, 
and the Sobolev inequality 

Y < (2/DN )(4/D - l)D/2-l X D/2 N 2-D/2 (7) 

This expression proceeds from optimizing the estimate 
Y < CX D / 2 N 2 ~ D / 2 by means of variational methods, 
which amount to minimizing the functional J{u} — 
X D l 2 ~N 2 ~ u I 2 jY and provide a dependence between the 
best constant Cb es t and the quantity Nq. Here, Nq is 
the value of N computed with the radially-symmetric 
ground-state Xo(f) satisfying the NLS equation — xo + 
V 2 xo + Xo = with no trap, such that No = 11.68 for 
D = 2 and iV = 18.94 for D = 3 @. When fl = [free 
NLS], a sufficient condition for collapse is H < 0. In the 
2D case, Eq. (0) then imposes that N must be larger than 
N CI = Nq ~ 11.7. In the 3D case, sharper estimations 
of the complete identity (§), d 2 (r 2 ) = (4/N)(3H - X), 
yield the more stringent sufficient condition for collapse 
H < Nq/N, after bounding H from below with (|]) 

When Slo ^ 0, the mean-square radius (r 2 ) satisfies 



where the strict equality concerns the 2D case only, and 
H < is also sufficient for collapse, since it leads to 
dti 1 " 2 ) < from Eq. (|^). Because H expands as 
H = Hq =q + H cr , where H cr = jNQ 2 R 2 , this require- 
ment implies Hq 0= o < 0. For D = 2, N must thus 
necessarily exceed N CI — Nq [plJl^]. In the opposite case 
H > 0, two distinct dynamics can develop: (a) When 
Hn =o < 0j collapse occurs again within the domain 
H < H CI for which the existence of a finite maximal 
blow-up time, tf 3 * = (l/2fi ) cos" 1 [1/(1- 2H CI /H)], 
makes sense. In the special case N = N cr (H = H CT ), 
2D wavefunctions blow up with an amplitude diverging at 
t c = ir/2n . (b) When ffn =Q > with N < N CT , there is 
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no collapse for D — 2: condensates having H > H CI never 
blow up and evolve by oscillating with the constant fre- 
quency 2f2 [@- For D = 3, the latter conclusion must, 
however, be subdued, because sharper criteria for blow- 
up may still apply ||,[l!|, so that H > H CI just consists 
of a necessary condition for the absence of collapse. 

Hence, although reinforced with a parabolic trap jl3) , 
the collapse generally manifests under specific conditions, 
such as H < H CI = jNHqRq. In the opposite case, 
H > _ff cr , collapse can disappear and the wavefunction 
u does not systematically spread out, but may instead 
oscillate or even form robust solitonlike states. This ■pos- 
sibility of forming stable condensates remains an open 
problem for high dimension numbers. Therefore, it is 
worth knowing whether the stationary bound states of 
Eq. (0), expressing as 



u(r,t) = x{r)exp(iAt), 



(9) 



can exist for spatial dimensions D > 2. These ground 
states are defined by the solutions of the differential equa- 
tion 



n 2 

-Ax + V 2 x-^r 2 X +\x\ 2 X = 0, 



(10) 



which functionally depend on A and are localized in 
space, i.e., x{^) = at r -> +oo. Here, the parame- 
ter A is sign-opposite to the energy eigenvalue E = —A, 
which corresponds to the above-defined chemical poten- 
tial /j. in reduced units. To start with, we examine the 
precise conditions under which these states may form. 
We multiply Eq. ( |To| ) by x* an( i ^ • ^X*> then combine 
the results to find 



X s = 
Y x = 



J ^- )[ DA + %D + A)(r 2 ) ! 



(11) 



(A-D) 



[2A + 2 (r 2 



where subscript s applies to the integrals computed on 
the ground-state solution. As X s and Y s are both posi- 
tive quantities, it is easy to check that localized solutions 
X exist for D < 4, provided that 



O 2 

A > ~f (r 2 ). 



(12) 



Unlike their free NLS counterparts (Sl = 0), confined 
stationary states can thus exist for negative parameters 
A and the requirement ( |T^ ) implies that the energy E 
must be less than twice the mean-square radius of the 
stationary condensates. Furthermore, the integral H for 
the ground states reads 



d - 2 , n 2 , 9l 

H s = AN S + ?—N s (r 2 )< 



(13) 



which ensures d 2 (r 2 ) s = from Eq. (||), as expected 



Combining relation ( |l3| ) with the existence condition (|1 
then yields the constraint that the Hamiltonian must sat 
isfy on the stationary solutions: 



H M > 



n 2 

fN s (r 2 ) 



> 0. 



(14) 



where 



Thus, H s for steady solutions of trapped BECs is strictly 
positive. Moreover, we can observe from the result (|lT 
that H s belongs to the range 
there is no collapse for D = 2 if N s < N cr , and where 
collapse can be prevented for D = 3. On the other hand, 
since collapse implies the blow-up of the gradient norm 
X — » +oo as (r 2 ) — > [see Eq.(|^)], a sharper condition 
assuring the collapse at dimensions D > 2 can be in- 
ferred in the domain X > X s from the virial relation (|^) 
rewritten as 

d 2 (r 2 ) = ^-[H-(l- 2/D)X] - (D + 2)^<r 2 ), 

and it reads H < (1 — 2/D)X s . From expression (|l3|), 
H s can be seen to never belong to this class of Hamilto- 
nians. Therefore, as H s always lies above this bound to- 
gether with the limit value H CI = ^Ntt 2 ,]! 2 , below which 
any initial datum certainly promotes the collapse, it is 
important to investigate the stability of the stationary 
states x- 



III. STABILITY CRITERION 



We first use the heuristic argument, following which 
the functional dependences of H on u under the con- 
straint of a fixed particle number N provides some qual- 
itative information about stability. Indeed, as is readily 
seen from the variational problem 5(H + AN) = 0, x 
governed by (^) realizes an extremum of H at fixed N. 
Under this constraint, we employ the compatible sub- 
stitution u = a~ D / 2 x{r /a), where a plays the role of a 
Lagrange multiplier, and we plug it into H to obtain 



H a = X s /a 2 



Y s /a D +a 2 I 3 , 



(15) 



where the integrals X s and Y s are here expressed in terms 
of £ = r/a, while I 3 = j^ 2 f £ 2 \ X (£)\ 2 dC K is then 
clear from Eq. ( |l5| ) that H a admits a global minimum 
for D — 2, whenever X s > Y s . Also for D = 3, H a 
may exhibit a local minimum under the same condition. 
The minima of H a are given by the roots of the iden- 
tity SH a /Sa\ a=1 = 0, i.e., X s - DY s /2 - 7 3 = 0, which is 
nothing else but the characteristic relation for the ground 
states, as this can be refound from Eq. (|Ti"l). This ar- 
gument points out that the minima of H are reached on 
the stationary states %, that constitute stable equilib- 
rium solutions around which any nearby solution can be 
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trapped, provided that X s > Y s . Note that this inequal- 
ity is always satisfied by all bound states x when D > 2. 
For comparison, with the same arguments the free NLS 
ground states (f2o = 0) are found to be unstable and H a 
has no global minimum for these dimension numbers. 

To derive a more rigorous proof for the stability of 
stationary wavefunctions, we now perturb the latter by 
means of the following solution 



u(r, t) — [x(r) + v(r, t) + iw(r, t)]e 



iAt 



(16) 



where v,w <C x are rea l functions, which we assume to 
be radially-symmetric. Inserting ( |l6| ) into Eq. (^J), lin- 
earizing and decomposing the resulting equation into real 
and imaginary parts then yield 



Lqw = dtv, L\v = —dtw. 



where the differential operators 



(17) 



Lo = -V 2 +A- X 2 + ^r 2 , L 1 =L Q -2 X 2 , (18) 



have the properties 

L 1 d X /dA=- x , (19) 
L t Vx = -ngfJc/2, L fx = -2V X . (20) 



vectors Vx and — r*x are both orthogonal to the ground 
state x- A combination of Eq. (EG) then leads to 



Xo^iVx = «qVx 



(22) 



and, thus, Vx is an eigenstate of LqLi < L\ with positive 
eigenvalue. For Lq being positive definite, this property 
indicates that Ai > 0. The lowest eigenvalue Ao of L\ 
is necessarily negative, although larger than for = 0, 
because < r\\L\T] > attains negative values with, e.g., 
t] = x- Moreover, L\ has exactly one negative eigenvalue 
in its spectrum [0. These technical conditions being 
fulfilled, we look for perturbations growing as v,w ~ e rt 
with growth rate T and determine the condition assuring 
the stability of ground states by maximizing 

<v\L lV > (23) 



r 2 = 



< v\L n L v > 



in the subspace of functions v orthogonal to x- Stability 
then results from the proof that < v\Liv > is positive 
under the constraint < v\ X >= 0. According to stan- 
dard procedures [jl8|-po"|l , we identify the discrete spec- 
trum of Li by setting L\v — \v + ax, where a is a La- 
grange multiplier related to the orthogonality condition 
< v\x >= 0. We construct a complete basis of orthonor- 
malized eigenfunctions as \v >= J2n Cn^n >, such that 
Li\ip n >= \n\ipn >, and combine these relations to get 



G(A) = ^ <X| ^ >< ^ |X> =I). 



A n — A 



(24) 



with Vx = (r/r)drX- It is interesting to notice that the 
small-amplitude perturbations described by Eqs. jl^ ) 
to (JlJ) provide an alternate form of the usual Bogoli- 
ubov excitations. Here, x is supposed to be the unique, 
radially-symmetric solitonlike bound state of Eq. JIc| ) 
with no zeroes. So, as L = — x _1 V[x 2 V(l/x)] is non- 
negative and admits the discrete eigenvalue for x only, 
zero is the lowest eigenvalue of Lq. In addition, since x 
has no node, Vx possesses a single node. We then apply 
a basic theorem from spectral theory claiming that if ipk 
is an eigenstate for L± < Lq such that L±ipi. — Xkipk, and 
if ipk has exactly k zeroes, then Afe is the (k + l)th eigen- 
value ranked as Ao < Ai < ... < Afe. This theorem holds 
when the eigenfunctions of L\ are L 2 -integrable, which 
is here ensured since the linear eigenstates of L\ decay 
much faster with Gaussian tails than their free NLS coun- 
terparts. For Qq = 0, Ai is equal to zero and, therefore, a 
unique negative eigenvalue Ao < exists. Let us now ob- 
serve what happens with a trapping potential (Qq 7^ 0). 
First, L\ has certainly a strictly positive eigenvalue since 



< Vxl^iVx >= 



> 0, 



(21) 



where < | > denotes the standard L 2 scalar product ap- 
plied to real functions. Second, it is obvious that the 



In (|4j), the function G(A) monotonously increases from 
— oo to +oo in the range ]Ao, Ai[. With Ao < and Ai > 
0, stability follows from the sign of G(0) =< xl-^i^X >> 
which must be negative in order to ensure A > 0. As 
G(0) is defined by G(0) = — < xl^x/^A >, we conclude 
that a sufficient condition for stability is given by 

This condition guarantees the orbital stability of the sta- 
tionary condensate. By "orbital" stability, it is meant 
that, modulo the elementary symmetries of Eq. (Q), the 
shape of a BEC soliton is preserved by perturbations hav- 
ing no growing modes, when they act on any solution 
staying nearby the ground state orbit fll7|] . Conversely, 
dN s /dA < is sufficient for any bound state of Eq.(||) to 
be unstable. If we moreover admit that dN s /dk — also 
leads to instability, condition is not only sufficient, 
but also necessary for the stability of ground states. 



IV. NUMERICAL RESULTS 

Identified through numerical integrations of Eq. (^) , 
some radially-symmetric, bell-shaped solutions x are pre- 
sented in Figure 1 for different parameters A and f^o = 2. 



4 



Two-dimensional and three-dimensional ground states 
are illustrated in Fig. 1(a) and Fig. 1(b), respectively. 
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FIG. 1. Examples of radially-symmetric, localized station- 
ary condensates with different A and Qo ~ 2, at both dimen- 
sions D = 2 (a) and D = 3 (b). 

The dependence N s versus A is shown in Figure 2 
for condensates defined at dimension numbers two [Fig. 
2(a)] and three [Fig. 2(b)], with fio = 2. Dotted lines in- 
dicate the variations of the ground-state particle number 
N 8 with respect to A for the free NLS equation (Sin, = 0), 
namely, N s = 11.68 in 2D and N s = 18.94/VA in 3D 
JL5| . Note that N s for fio / always lies below these 
free particle numbers. In the limit A — > +oo, the solu- 
tion x behaves as free NLS solitons, which are all unstable 
with dN s /dh < 0. In contrast, as can be seen from this 
set of figures, trapped ground states are stable for every 
A with moderate values when D = 2, whereas for D = 3, 
the criterion (|25| ) predicts their stability for negative A 
only, i.e., A < A cr , where A cr ~ —0.72 corresponds to the 
maximum number 7V] nax ~ 14.45. 




-5 5 10 -5 5 10 

A A 

FIG. 2. N s vs. A for the ground states of Eq. (2). Solid 
curves refer to 2D (a) and 3D (b) confined ground states. 
Dotted lines indicate their free NLS limits. 

The difference between "free" and trapped solutions 
is that the discrete spectrum of L\ — + fl^r 2 /4 

is here shifted towards the range of positive eigenvalues. 
Thereby, stability of ground states seems "reinforced" 
by the parabolic trap, compared with the case fio = 0. 
Such a result cannot be obtained for a detrapping po- 
tential (fig — > — fig), because the linear eigenstates of L\ 



are not L 2 -integrable in this configuration. Even if the 
above procedure could be applied in that case, Eq.([2l|) 
furthermore suggests that at least two negative eigenval- 
ues would appear in the discrete spectrum of L\, from 
which instability surely follows. 

Finally, we have plotted in Fig. 3 the temporal evolu- 
tions of solutions to Eq. (^) with a parabolic trap for var- 
ious initial data. In the 2D case, Fig. 3(a) shows |tt(0, t)| 2 
for initial conditions defined by the ground states x with 
A = 1 and A = 5 for fig = 2, and with A = 1 for 
fio = 0. In the 3D case, Fig. 3(b) represents similar evo- 
lutions from ground states having A = —2, —1,0, among 
which the solution initiated with A = inexorably col- 
lapses. Dots show the same initial data undergoing small- 
amplitude periodic perturbations. From those, we can 
observe the robustness of the 2D and 3D ground states 
of the GP equation, for which both inequalities ( |l4| ) and 
( p5| ) are always verified. Note that perturbations do not 
affect the ground state orbits, as long as A is far below the 
critical values beyond which stationary solutions must 
become unstable, i.e., A 3> 5 for D = 2 and A > —0.72 
for D — 3. Stability is lost as A attains large positive 
values, in accordance with the criterion (B5j) . 
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FIG. 3. Temporal evolution of |u(0, t)\ 2 for different initial 
data with N = N s . (a) D = 2, A = 1 (solid curve), A = 5 
(dashed) for Qo ~ 2, and A = 1 for fin, = (dash-dotted). 
This datum promotes a finite-time collapse at t c = n/2Qo- 
(b) D = 3, A = -2 (solid), A = -1 (dashed), and A = 
(dash-dotted) for Qq — 2. The dots illustrate the action of 
perturbations introduced at t — around the ground states. 



V. CONCLUSION 



We have discriminated the different regions of 
collapse/no-collapse for the non-stationary solutions of 
Eq. (^) and determined the exact conditions under which 
stationary localized states of this equation can exist. We 
have also demonstrated that the stationary ground states 
lie in the region H > H cr > 0, where collapse is defini- 
tively absent in 2D and can be avoided in 3D, and that 
they are stable provided dN s /dA is positive. This re- 
sult displays evidence that a parabolic trap makes the 
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multi-dimensional, stationary solitary modes of the GP 
equation stable for a wide class of parameters accessible 
in BEC physics. 

To conclude, the criterion (|2^) allows for determining 
straightforwardly the stability of condensates from their 
particle numbers only, which should be useful in current 
experiments. 
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